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The established macroscopic equations of motion for two phase im- 
miscible displacement in porous media are known to be physically in- 
complete because they do not contain the surface tension and surface 
areas governing capillary phenomena. Therefore a more general system 
of macroscopic equations is derived here which incorporates the spa- 
tiotemporal variation of interfacial energies. These equations are based 
on the theory of mixtures in macroscopic continuum mechanics. They 
include wetting phenomena through surface tensions instead of the tra- 
ditional use of capillary pressure functions. Relative permeabilities can 
be identified in this approach which exhibit a complex dependence on 
the state variables. A capillary pressure function can be identified in 
equilibrium which shows the qualitative saturation dependence known 
from experiment. In addition the new equations allow to describe the 
spatiotemporal changes of residual saturations during immiscible dis- 
placement. 
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I. INTRODUCTION 



An aquifer or a petroleum reservoir may be seen as mixture of porous rocks and 
soil with various gases and liquids, usually residing at depths that prohibit precise 
measurements on the constituents 

UH- Despite many years of research and despite 
the fact that the microscopic laws of fluid dynamics are well known the prediction 
of macroscopic multiphase flow in porous media has not met with success (see 
for recent overviews). 

Microscopic and macroscopic descriptions of multiphase fluid flow in porous me- 
dia differ considerably from each other and both have their own characteristic prob- 
lems Jl]J^||]. A microscopic description fails because it is generally impossible and 
not interesting to know the detailed microstructure and flow patterns on the pore 
scale. In contrast herewith the generally accepted macroscopic models are incom- 
plete. Owing to its widespread importance in applied problems the transition from 
the microscale to the macroscale for two phase flow in porous media has attracted 
considerable interest from physicists in recent years P-[2^|. Relatively little or no 
attention has been focussed, however, on the continuum limit of the proposed mi- 
croscopic models. Effective macroscopic equations of motion used by engineers are 
based on generalizations of Darcy's law, although every practitioner knows that they 
are incomplete [^3|,|fJ . Microscopic models by physicists often attempt to predict the 
relative permeabilities and the capillary pressure functions for simple microstruc- 
tures although it is clear that they are highly nonunique. 

Detailed and rigorous derivations of macroscopic equations from microscopic ones 
is usually the last step in understanding a physical phenomenon. Experimental 
observations and tests of phenomcnological macroscopic models usually precede 
the microscopic understanding. It is therefore of interest to improve the current 
macroscopic description which is known to fail even for simple laboratory systems. 




Given the need for more predictive macroscopic theories this paper proposes a 
new set of macroscopic equations. Loosely speaking the traditional theory includes 
only mass and momentum balance but fails to incorporate the energy balance. 
Obviously the interfacial energy depends on the interfacial area and hence the set 
of macroscopic observables describing the state of the mixture must be enlarged 
to include surface areas. Relative permeabilities and capillary pressures must then 
be identified within the new framework. I shall begin the discussion by presenting 
briefly the traditional formulation and discussing its problems. Afterwards the 
macroscopic balance laws for mass, momentum and energy will be combined with 
constraints and constitutive relations to give new macroscopic equations of motion. 
My last section in this paper will discuss the consequences and the identification of 
relative permeabilities and capillary pressures. 



II. PROBLEMS WITH EXISTING MACROSCOPIC EQUATIONS 



Consider simultaneous flow of two incompressible and immiscible fluid phases 
denoted generically as water (subscript W) and oil (subscript O) inside a porous 
medium. The accepted and widely used macroscopic equations of motion pU4|iap3] 
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are based on mass conservation 
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and the well known generalization of Darcy's law |l|,[?5| 
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all of which are assumed to be valid on length scales large compared to a typical pore 
diameter. The variables in these equations are the macroscopic pressure fields of 
water and oil, denoted as P w and P a , the water and oil saturations S w , S a and the 
macroscopic (Darcy) velocity fields v w , v . The saturation is defined as the ratio 
of water (oil) volume to pore space volume. Pressures and saturations are averages 
over a macroscopic region much larger than the pore size, but much smaller than the 
system size. Their arguments are the macroscopic space and time variables (x, t) . 
The porosity (f> is the ratio of pore volume to total volume of the medium. The fluid 
viscosities and densities are denoted by p, and p. The acceleration of gravity is g and 
the function Z(x) is the depth function. In the equations above K stands for the 
absolute (single phase flow) permeability tensor, k r w is the relative permeability for 
water, k r a the oil relative permeability. The equations above are closed by defining 
the capillary pressure as 



P c (x,f) =P (x,f)-P w (x,t) 



(4) 



and by postulating for it the constitutive relationship 



P c (x,i) =P c (S n (x,i)) 



(5) 



where 



S n (x,t) = 



>Sw (X, Sy 

1 S w i — So 



(6) 



is a normalized saturation. The water saturation obeys S W i < S w < 1 — 5 or 
where the two numbers < S W i , S or < 1 are two parameters representing the 
irreducible water saturation, S W i, and the residual oil saturation, S or . The residual 
oil saturation gives the amount of oil remaining in a porous medium after a water 
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flood. The normalized saturation S n varies between and 1 as S w varies between 
Syyi and 1 $ ov • 

Next it is argued [^,0 that each fluid flows in flow channels given by the solid 
matrix and the presence of the other fluid. One postulates that to each saturation 
there corresponds a unique configuration of flow channels. This picture is then 
formalized into the constitutive assumption that the relative permeabilities 



k r w (x,t) = k r w (S n (x,t)) (7) 
fc£(x,f) = *£(S„(x,t)) (8) 

are functions of saturation alone. 

Although these equations are widely accepted and almost universally applied in 
numerical reservoir simulation ]27j ] or aquifer modeling [p8| they must be considered 
physically incomplete. Even if the picture of flow channels is accepted it is clear that 
the resulting permeability depends not only on saturation (i.e. volume fraction) but 
also on the surface area with other fluid or solid phases. This experimental fact is 
well known from observations of single phase flow in porous media [p6jj 



The second basic problem with the macroscopic equations of motion (|])-(H) arises 
from the experimental observation that the parameters S or and S W i are neither 
constant nor are they known in advance for a given experiment. Instead they 
vary with space and time and depend on the flow conditions in the experiment 
|29|-|37|l. More precisely, the residual oil saturation depends on the microscopic 
capillary number Ca = [J. w v/a ow where v is a typical flow velocity and <7 ow is 
the surface tension between the two fluids. The dependence S or (Ca) is called a 
capillary desaturation curve. The measured values of S or are known to depend also 
on the duration of the flood. Capillary desaturation experiments contradict clearly 
to the assumption that the functions k^,(S w ),k^,(S w ) and P C (S W ) depend only 
on saturation. Instead they depend also on velocity (3^]3!| and pressure. Hence 
they depend on the solution and cannot be considered to be constitutive relations 
characterizing the system. Therefore the system of equations of motion must be 
considered to be incomplete. 



III. MACROSCOPIC CONSERVATION LAWS 



The basic idea of this paper is to include the energy balance involving the macro- 
scopic interfacial energy into the basic equations of motion. As the interfacial 
tensions do not appear explicitly in equations (0) through (0) it is clear that this 
requires an enlarged description of the macroscopic state. While previously the 
macroscopic state of the fluids was described by saturations, pressures and veloci- 
ties it must now include macroscopic interfacial areas per unit volume. The necessity 
to include interfacial areas in a macroscopic description is also stressed by volume 
averaging approaches [^0]-^2j . The main difference between volume averaging tech- 
niques and the approach of the present paper is that the present approach is based 
purely on macroscopic mixture theory, and makes no reference to the underlying 
microscopic equations. 
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To discuss mass, energy and momentum balance laws in a unified manner consider 
first a general conserved quantity W defined per unit mass. Let r denote the flux of 
this conserved quantity across a surface of a material control region &(t) in R 3 . The 
general multiphase conservation law in the theory of mixtures reads in differential 
form plEl 



— (^p.fj + 4> a p a ^ a V • v a - V • t„ - 4> a p a g a = e Q (9) 



Here the subscript a refers to the different phases of the mixture, and 

D a /Dt = d/dt + v a -V (10) 

denotes the material derivative and v Q is the velocity. p a is the mass density, and 
4> a is the volume fraction of phase a. Of course 

and introducing the volume fraction of the void space, i.e. the porosity, 

a^solid 

allows to define saturations of the fluid phases through 

0„ - 05„ (13) 

In (^|) (7 Q denotes the external supply of ^ . e a represents the interphase transfer of 
the conserved quantity VP into phase a from all other phases. In other words 



where e Q/3 is the transfer of \& from phase j3 into phase a. Of course global con- 
servation of 4* would require that summation of (j^) gives zero, e a = 0, but 
this will not be required here. This general formalism can now be applied to mass, 
momentum and energy conservation. 

For mass conservation f„ = 1,t q = 0, g a — and e Q = M a is the transfer of 
mass from all other phases into phase a. Inserting into eq. (|9|) yields 

^(<t> aPa ) + V -{cj> aPa v a )=M a (15) 
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as usual. 



For momentum conservation — v Q is the velocity, t q = S„ is the stress tensor, 
g a = b Q is the body force, and e„ = m Q is the momentum transfer from all other 
phases into phase a. Inserting into (^J) and using mass balance now yields 

(f> a p a - V • S„ - (f> a p a b a =m a - w a M a (16) 

where the divergence V • S a of the stress tensor is a vector whose i-th. coordinate 
is Ytjdi^ij/dxj. 

For energy conservation only the kinetic energy and the interfacial energies are 
considered, while thermal effects will be ignored. Hence in this case 



2 2^^+^p, V ; 



where a a/3 is the surface tension between phase a and (3 and A ali is the specific 
internal surface of the a/3-interface. The denominator in the interfacial energy 
terms (giving the a/3-density) ensures the same normalization with respect to mass 
and volume as for the kinetic energy term. The factor 1/2 in front of it attributes the 
interfacial energy between phase a and p equally to the two phases. This assumption 
will be relaxed in a future publication pa]. Proceeding with this assumption the 
energy flux becomes 



(18) 



where S Q • v„ is the work done by stress and q Q is the flux of surface energy. The 
external energy supply is similarly split, 



(19) 



into the work b Q • v Q done by body forces and the external supply h a of interfacial 
energy. Finally the energy exchange from all phases into phase a is e Q = E a . 



Inserting equations (|17|) through (|19| ) into the general balance law (|9j) gives the 
energy balance for immiscible displacement as 



Vv„ — 4>c,p a h a 



= E a -m a -v a + M a -£ (20) 
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TABLE I. Overview over the subphases and their indices for immiscible displacement 
of two fluids inside a rigid porous medium. 



Phase index a 


Phase description 


1 


connected (percolating) water phase 


2 


disconnected (trapped) water phase 


3 


connected (percolating) oil phase 


4 


disconnected (trapped) oil phase 


n 


rigid porous rock phase 



where 

d d p., \ 

8=1 j = l 

Equations (|l7|) and (^(j|) embody the idea to include interfacial areas and energies 
into the macroscopic energy balance. This is the first step beyond the traditional 
approach described e.g. in p] ]26|j43ll . 

The second step beyond the traditional theory consists in splitting each fluid 
phase into a connected (percolating) and disconnected (nonpercolating) subphase. 
The idea here is that such a division is necessary to handle the dynamics of residual 
phase saturations. During an immiscible displacement process a finite fraction of 
each fluid phase becomes trapped in the pores. Once trapped the droplets are 
immobile and fixed through large local capillary forces. Only at very high flow 
rates can such trapped droplets be mobilized. Therefore the second main idea of 
the present theory is to treat these droplets as an individual phase interacting with 
the connected mother phase mainly via mass exchange. 

Consider now the simplest case of two immiscible fluids (called oil and water 
as before) flowing inside a rigid porous medium. Dividing each of the two fluid 
phases into two subphases one has a total of five phases. Each fluid superphase 
is divided into a percolating (connected) and a trapped (disconnected) subphase. 
Here percolating means that each point in the subset occupied by the percolating 
subphase can be connected to the exterior boundaries of the sample by a path 
within the subphase. The disconnected subphase is the complement within the 
chosen superphase. The connected (percolating) subphase can also be defined as 
the region inside of which an external applied pressure gradient can propagate. 
Therefore the phase index a above can assume five values summarized in Table Q. 



IV. CONSTITUTIVE RELATIONS AND CONSTRAINTS 



This section complements the basic balance laws with geometrical and physical 
constraints and with constitutive relations. 
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The fluids are assumed to be incompressible 



Pl(x,t) = p 2 (x, t) = p w — const (22) 
p 3 (x, t) = p 4 (x, t) = /9 = const (23) 

and their stress tensors are approximated as 

S Q = -0 Q P Q 1 (24) 

for all a G {1, 2, 3, 4} where 1 is the unit tensor. In ( pi| ) it is assumed that all shear 
stresses in the fluid are negligible, and that P a are the true pressures in the fluids. 
It will be assumed that the only body force is gravity, b Q = g VZ where g is the 
acceleration of gravity and Z(x) is the depth function. In addition inertial effects 
will be neglected, D a v a /Dt = 0, for all phases. All of these assumptions are also 
made in the derivation of the traditional equations described in the previous section 



It will be assumed that there are no external sources or sinks of surface energy, 
g a = 0. Next it is assumed that q„ = 0, i.e. there is no diffusive flux of surface 
energy during the time scale of observation. A crucial assumption is that the ex- 
change of interfacial energy between the phases is dominated by separation and 
coalescence of disconnected subphase droplets from/with their connected mother 
phase, i.e. through mass transfer, and that all other mechanisms are negligible. 
This amounts to setting E = 0. 

The porous matrix as well as the disconnected (trapped) subphases are immobile, 
and hence 

v K = v 2 = v 4 = 0. (25) 

The interfacial area between a disconnected subphase and its superphase vanishes, 
and the same holds for the contact between the two disconnected subphases. Hence 

A12 = A 34 = A 24 = (26) 

Of course A aa = and symmetry A afj — Ap a holds. 

The disconnected fluid phases have macroscopic pressures P 2 (x, t) and P 4 (x, t) 
obtained by averaging the pressure in the disconnected drops and droplets. The 
resulting functions -P 2 (x, t), Pa(x., t) however are generally not continuous and hence 
not differentiable. 

The wetting properties are assumed to be constant, and they are specified by the 
requirements that the individual interface areas with the rock are constant, i.e. 

(27) 
(28) 
(29) 
(30) 



Aj-ji (x, t) 


= Am 


= const 


A 2 *(x,t) 


= A 2n 


= const 


A 3K (x,i) 


= A 3n 


= const 


A 4TC (x,t) 


— a 4tl 


= const 



7 



Of course A = A\ n + A 2 -^ + v4 3k + A^ is the total specific internal surface of the 
porous medium. The ratio 

Ww = A ^ + A ^ (31) 

varies from 1 to as the medium varies from water wet to oil wet. 

Finally it remains to specify the momentum exchanges m,-, = ^2 m a/3 and the 
mass exchanges M a — ^2 /3 M a/3 . In both cases m QQ = 0, M aa — will be used. 
The momentum exchange is conventionally modelled [ fl3"| through Stokes drag in 
the form 



nw = R a0 (v - v„) = — -(v„ - v a ) (32) 

where the mobility A Q|3 = fc Q(3 //z Q/3 is expressed as usual in terms of the permeability 
k aft and the viscosity /j, al3 . As a new assumption the mobilities A Q|3 are modeled 
through 

A - = ¥ = S- (33) 

Ma/3 ^0^0/3 

analogous to the Carman-Kozeny relation for absolute permeabilities. This assump- 
tion reflects the traditional view that each fluid flows in its own channel. If there is 
no common interface, then it is not possible to exchange momentum, hence 

mi2 = m 34 = m 2 4 = ( 34) 

completes the specification of the momentum exchange. 

The rigid porous matrix has a volume fraction 1 — =const which is assumed 
to be spatially and temporally constant. The rigid matrix is dynamically inert (see 
eq.(p5|)), and hence only the balance laws for the fluid phases will appear below. 1 

Chemical reactions between oil, water and rock are excluded, and hence 

Mi 3 = M 14 = M 23 = M 24 = M ln = M 2n = M 3K = M 4n = 0. (35) 

Mass transfer occurs within each fluid phase between its connected and the discon- 
nected subphases. These nonvanishing mass transfers are 

M W = M 12 = -Mil (36) 



1 Because the rock matrix is assumed to be rigid and impermeable one has Ftjn = R-n2 = 
R-r.3 = Rni = oo which implies dissipation. Because of this and (^) the balance laws for 
the matrix phase would either be trivial or ill defined. 
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M = M 34 = -M43 (37) 



In general the mass transfer may be a function of the velocities, saturations 
and interface areas of the other phases. A precise form for the dependence 
M w = M w (vi, A13, Si S3) is not needed here. Any concrete assumption for such a 
dependence would have to be checked by experiment. 



V. RESULTS 



Combining the balance laws with the constitutive equations gives the macroscopic 
equations of motion of the two pore fluids. 2 Mass conservation yields four equations 

d 

Pw <j>—S x + Pw (f>V • (Sivi) = M w (38) 

p w ^S 2 = -M w (39) 

Po ^S 3 + Po <t>V -{S 3 w 3 ) = M (40) 

Po ^Sa = -M (41) 

while momentum conservation gives only two equations, 



V{4> 1 P 1 )-^ lPw gVZ = RisVa-iRis + Ru + Rm + M^V! (42) 
V((/> 3 P 3 ) - foPog VZ = i? 31 vi - (R 31 + R 32 + Rsn + M )v 3 (43) 



because P 2 and P4 are not differentiable. Energy conservation gives again four 
equations 



1 + ( 4>lPw<7 wo A 13 ^ 4>lPw<T wo A 1A ^ 01/O w cr WK A l7 



2\ Dt J \Pw4>i + Po03 Pw4>i + Po4>i Pw0i + Pk(i-<; 

v 2 

= -<AiPi V • vi - P13V1 • v 3 + (i?i 3 + RiA + i?i K )v 2 + M w ^- (44) 



9 ( (j)2Pw^woA 23 , 02 Pw f7 'w7^^2T 



dt \ p w (j) 2 + Po4>3 Pw4>2 + Ptz(1 



= (45) 



2\ Dt J \p m 4>i + p 4> 3 p w (j) 2 + Po4>3 Po^3 + Pw(l 



2 Balance equations for the rock matrix do not arise (see footnote 1). 
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3 P 3 V ■ v 3 - i?3ivi • v 3 + (i? 3 i + R32 + R3k)v 2 3 + M y (46) 



8 ( (t>4Po<T wo A u 04 



dt \p w (t>l + po4>4 Po4>4 + Pn(l - ■ 



= (47) 



The unknown fields in these equations are the four saturations S a (a = 1,2,3,4) 
of the fluids related to </> Q by equation (|l~3|), the two pressures Pi, P3, the six com- 
ponents of the two velocities vi, v 3 of the connected fluids, and the three surface 
areas A13, An, A 2 3- Note that all A a ^,a — 1,2,3,4, are constant by virtue of 
eqs. d27j ) through fl30|). This is a total of 15 unknowns for which the condition 
Si + S 2 + S 3 + 5 4 = 1 (see Q)) together with equations (||)-© form a sct of 15 
equations. 

A detailed analysis of these equations will be presented elsewhere 0] . Here only 
the simplest comparison with the traditional equations (jl])-(||) will be discussed. In 
particular the question arises whether the new equations allow to determine relative 
permeabilities and capillary pressure functions. 

The capillary pressure, defined in eq. (|4|), is an equilibrium property. In the 
traditional equations it is related to the saturation profile in complete gravitational 
equilibrium p6| , i.e. when v„ = and when all terms involving time derivatives 
are set to zero. Proceeding in this way one finds that the resulting equations can 
be solved to give 

P c (x,oc) =P 3 -Pi= C3 - — Cl + (po - Pw)x • g (48) 
0D 3 (X, 00) (pi>i(x, 00) 

where g is the acceleration vector of gravity and c 3 ,ci are integration constants. 
Using S w = S1+S2 and 1 — S w = S3 + S4 shows that in the directions perpendicular 
to gravity 

P c (x, 00) = P 3 - Pi = 03 - Cl (49) 

<p{l-b4-b w ) (p{b w - b 2 ) 

If it is assumed that S 2 and Si are spatially constant throughout the sample then 
it is tempting to identify them with S W i and S^,. because this would imply the as- 
sumption (Q) of the traditional formulation. Moreover the capillary pressure would 
indeed have the familiar S-shaped form as a function of saturation. It should be no- 
ticed however that eq. ( fi"9| ) holds only without flow while assumption (||) is assumed 
to hold always. Therefore in the new formulation the capillary pressure depends 
strongly on the flow regime. This prediction agrees with what is known from theory 



and experiment | p8|]39|l . Notice also that in (49) the quantities 5*2 and S4 are state 
variables and hence the capillary pressure includes the spatiotemporal variability of 
the residual saturations. 

The generalization of Darcy's law follows from the equations of momentum bal- 



ance (W2f) and (43) by solving these equations for vi and V3. This yields a generalized 



form of Darcy's law for two phase flow in the form 

k k 
v i = fc ww — (V(</>iPi)-</>ipwg)+A^ — (V(4> 3 P3)-hpoS) (50) 
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k k 

v 3 = k r ow (yfap L )-fa Pw g) + kl — (V(0 3 P 3 ) - ^3Pog) (51) 

Mw Mo 

where the relative permabilities k r a „ are defined by these equations. The k r ryH are 
related to the resitivities R a/3 which in turn are given by equations ([52|) and (p3|). 

This result differs from the traditional gen eralization of Darcy's law given in eqs. 
(0) in three respects. Firstly ( ^p|) and ((5^) includes off-diagonal coupling terms 
k 7 wo and k r ow which are missing in (|). Secondly (|c]) and (|l]) include a new 
driving force 0P Q V<S Q resulting from the term V(</> Q P Q ). The new driving force is 
proportional to saturation gradients. Thirdly, and perhaps most importantly, the 



relative permeabilities defined through (|50| ) and (51) depend not only on saturation 
as assumed in (||) but 

Kf, = C(- R 13^31,i?14,i?23,i?l K ,i?3^,M w ,M ) (52) 

= k r JSi,S 2 ,S 3 ,S i ,A 13 ,Au,A 3 2,dS 2 /dt,dS i /dt,A 11 i,A3Ti) (53) 



by virtue of (t*2j) , (|43|) and (g2J),(|3J),(g9|) and (|4l|). This expresses a dependence on 
saturation, changes of saturation, disconnected (residual) saturations, interfacial 
area and wetting properties (Ai-^, A 3 ^). Thus the present theory gives a much 
more detailed picture of relative permeabilities than the traditional assumption (|J) . 

In conclusion the present paper has discussed the problems associated with the 
macroscopic equations of motion for two phase immiscible displacement in porous 
media. After discussion of the traditional formulation a new theory has been devel- 
oped based on mixture theory. The new equations differ from the traditional ones 
in the following main aspects: 



1. The set of macroscopic state variables is enlarged to include specific surface 
area of the fluid phases. 

2. The set of macroscopic balance laws is enlarged to include interfacial energy 
balance. 

3. The theory allows a macroscopic quantification of the wetting properties 
through the surface areas. 

4. The dynamics of the disconnected fluid phases is included. This allows to 
describe the spatiotemporal behaviour of residual saturations. 

5. The theory predicts a dynamic capillary pressure which, at least in complete 
gravitational equilibrium, shows the correct qualitative behaviour as a func- 
tion of saturation. 

6. The relative permeabilities are found to depend on macroscopic state variables 
other than saturation. In particular they depend on wetting properties and 
the internal surfaces areas. 



All of these results are consistent with experiment. More analysis of the equations 
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is necessary to compare their solutions with experiment. 
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